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1. Introduction 

Nonlinear Hamiltonian lattices are systems which are easy to formulate and to simulate 
numerically, but nevertheless their dynamics is far from being fully understood. 
Generally such a lattice is a non-integrable high- dimensional system demonstrating 
chaos. So one expects an effective "thermalization" to a statistical equilibrium. 
However, already the first attempts to follow such a thermalization initiated by Fermi, 
Pasta, and Ulam revealed many extremely nontrivial effects, still not completely 
understood (see Refs. [TJ [2] for recent progress of the FPU problem). While in the 
FPU model one is interested in the relaxation to equilibrium from an initial state which 
is localized in the space of modes (typically one or several modes in a finite lattice are 
initiated and one follows the transition to an equipartition), recently a problem of a 
spatial spreading of an initially local in space field has attracted a lot of attention in the 
studies of nonlinear lattices. Especially nontrivial is such a spreading in the case when 
linear propagating modes (phonons) do not exist. Then the spreading is due to nonlinear 
effects only. One way to suppress phonons is to consider lattices with disorder, in this 
case, due to Anderson localization, linear eigenmodes are exponentially localized and the 
spectrum is purely discrete. Recent numerical experiments with nonlinear disordered 
lattices have demonstrated that the initially localized wave packets spread in a very 
weak, subdiffusive manner [31 IH El El El El El HQl [HJ [121 H31 HH UHl H6l HZl HE]. The 
nontriviality of this observation is that the spreading is due to a weak chaos in the 
lattice, and it is not clear if this chaos survives when the density of the field becomes 
very small. All attempts to find the asymptotic properties of this spreading at very 
large times have been not successful so far, due to numerical problems of following the 
solutions for very long times. 

In this paper we develop the scaling approach to this problem, first suggested in [16] . 
We will try to establish and to check numerically the scaling properties of spreading, 
in dependence on the total energy of the initial wave packet. Our guiding tool will 
be the nonlinear diffusion equation which yields a subdiffusion spreading of a field. 
We formulate the scaling relations based on this equation, and check their validity for 
nonlinear lattices. 

We will start with formulating the object of our study, strongly nonlinear 
Hamiltonian lattices. We present a phenomenology of energy spreading and define the 
statistical quantities characterizing it (Section [2]). Next, we introduce a phenomenologic 
model that we use to describe the properties of the spreading process, namely the 
Nonlinear Diffusion Equation (NDE) (Section |3]). From its self-similar solutions we 
derive the scaling relations. In Section H] we present extensive numerical calculations for 
different types of nonlinearity. These results are compared with the predictions from 
the NDE. We end with some concluding remarks and outlook. 
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2. Strongly nonlinear lattices 

The main goal of this paper is to study properties of energy spreading in strongly 
nonlinear lattices. By strongly nonlinear we understand lattices where the coupling is 
described by nonlinear functions that disappear in the linear limit. So there are no linear 
waves (phonons) in such lattices and all propagating modes are nonlinear ones. Such 
lattices can be introduced in the framework of equations for the complex amplitudes 
(and then one obtains a strongly nonlinear generalization of the nonlinear Schrodinger 
equation) or as a generalization of a Hamiltonian Klein-Gordon lattice. In this work we 
follow the latter way. Moreover, we restrict ourselves to pure power-law nonlinearities. 



2.1. Hamiltonian 

In one dimension we formulated a strongly nonlinear lattice in terms of a Hamilton 
function for positions q k and momenta p k of oscillators labeled by site index k: 

H = J2^ + W^ k + % k+1 -q k ) X . (1) 

k K 

Here k>2 and A > 2 denote the powers of the on-site potential and the coupling term, 
respectively. For k = 2 we have a chain of nonlinearly coupled linear on-site oscillators; 
for k > 2 the on-site oscillators are nonlinear as well. Below we study situations with 
and without disorder, the latter is introduced via the variations of the parameters of the 
local potential u k (these are linear frequencies of the oscillators if k — 2 and parameters 
of the nonlinear on-site potential if k > 2). 



2.2. Reseating 

Hamiltonian (CQ) contains two parameters W and a that determine the time scale and 
the ratio of local to coupling potentials. For different powers of local and coupling 
nonlinearities, i.e. for k ^ A, we can get rid of these two parameters by rescaling the 
canonical variables and time as follows: 

q k -+ W b a- b q k , p k -> W Xh ' 2 a~^ 2 p k , t -+ W ^ b ' 2 a^ b 'H , (2) 

with b = — k). W and a then disappear from the equations and we are left 
with the total energy £ as the only relevant parameter depending on the initial state. 
Additionally, a distribution of local "frequencies" u k is relevant, while the width of this 
distribution is rescaled together with W. We will consider three cases: (i) no disorder, 
u k = 1; (ii) "soft" local oscillators, in this case the "frequencies" u k are chosen iid. from 
[0, 1]; (iii) "hard" local disorder, here the "frequencies" u k are chosen iid. from [0.5, 1.5]. 
In the rescaled coordinates the Hamiltonian reads 

Y2 T a • (3) 

k 

One special and highly interesting case occurs when the on-site and coupling terms 
have the same nonlinearity k = A. As now all terms containing coordinates q k have the 
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same power, one cannot set both parameters W and a to one by rescaling as before. 
Instead, one can use the remaining freedom to set the total energy £ to unity: 

q _> E 1/K W- xlK q, p^S 1/2 p, t W - 1/K £ 1/K ~ 1/2 . (4) 

Particularly, this means that the energy is not a relevant parameter of the system, but 
can rather be scaled to, say, 8 = 1, with a corresponding change of the time scale. Notice 
that the only remaining parameter is the ratio of strengths of on-site and coupling terms 
(3 = a/W. The rescaled Hamiltonian now reads 

H = ^Pl + M + ^-^\ (5 ) 
2. 3. Phenomenology of energy spreading 

For the Hamiltonian systems (15fo"j) we state the following question: How does an initially 
localized field spread over the lattice? We focus on very large systems, where boundary 
effects are not important. The distribution of the energy is characterized with its density 

*>k = Y = 8 ~ l (f + w t** + |f [ (gfc+1 ~ qk)X + {qk ~ gfc - l)A l) • (6) 

We start typically with non-zero values of Wk in a small interval (in most runs 10 sites), 
by choosing initial momenta randomly from a Gaussian distribution, and follow the 
distribution Wkif) in time. 

In the case of a lattice without disorder, regular waves can propagate along the 
lattice. Such localized solitary waves - compactons - have been thoroughly studied 
in [19] for a lattice with W — (i.e. without local potential). An initial localized 
perturbation emits compactons that dominate the process of energy spreading. For 
W 0, it is not known if exact compactons exist in a lattice. In numerics, we quite 
often observe "quasi-compactons" that propagate ballistically over large distances but 
lose energy and therefore eventually stop. We illustrate this in Fig. Ufa). In panel 
(b) we show the same initial conditions in a disordered lattice, here the propagation 
of "quasi-compactons" is blocked by disorder and one observes a slow spreading of the 
wave packet. This slow process can be described as follows. At each instance of time, a 
finite number of oscillators is excited while others have negligible energy. This effectively 
finite-dimensional Hamiltonian dynamical system demonstrates chaos. Due to this, the 
next neighbors to the excited region experience a weak chaotic forcing, and their energy 
slowly grows. After some time their energy reaches the level of the driving neighbor and 
they become involved into excited, chaotic domain, which in this way increases by one 
site. This site-by-site propagation means that in disordered, strongly nonlinear lattices 
at any finite time the distribution of energies is strongly localized, and has sharp edges 
(we expect that generally the field at the edges decays superexponentially fast, as for 
breather solutions in such lattices [20]). This sharpness is illustrated in Fig. [2j 
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Figure 1. Time evolution of an initially localized state for k — 4, A = 6, W — (3 — 1 
and energy £ = 10. (a): regular lattice, (b): lattice with "hard" disorder. The color 
coding corresponds to the logarithm of the local energy excitation log 10 Wk ■ The initial 
excitation was uniform on 10 sites. 
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Figure 2. Spreading of an initial single site excitation for k = 4, A = 6 and energy 
£ = 1.0. The plot shows local energy Wk vs. lattice site k for increasing times 10 4 , 10 6 , 
10 8 (inner to outer curves). Note the logarithmic scaling of Wk and the exponential 
drops in this scale. 



2.4- Measures of spreading 

In a statistical context, the spread of a distribution w k is measured via entropies, most 
suitable are the Renyi entropies 

1 y k 

that allow one to characterize also the spikeness/flatness of the distribution (in the 
context of energy spreading in disordered lattices this approach was introduced in [H]). 
We restrict here to the the entropies Ji and I2, which are nothing else than the usual 
Boltzmann entropy and the logarithm of the participation number P: 

h = -J2w k \nw k I 2 = -In^wl = InP . (7) 

k 

The participation number is a characteristic of the width of the wave packet rather 
popular in the context of Anderson localization studies [SJ El EI]. Both entropies 
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Figure 3. Schematic plot showing the two ways to measure spreading. Averaging at 
hxed L means averaging at fixed energy density, contrary to averaging at fixed time. 

define the effective width of the wave packet as = exp^ 2 ) (in particular, £ 2 — P)- 
As individual dependencies C(t) demonstrate enormous fluctuations, we perform an 
averaging of the entropies Iipit) over many realizations of disorder, thus obtaining 
smoothly growing widths C{t). 

For the strongly nonlinear lattices another approach is even superior to the 
calculation of the entropies. Here we determine the width L of the wave packet as the 
distance between its sharp edges as seen in Fig. [2] (independently on the distribution of 
the energy between these edges). After determining the spatial extend L, we measure the 
time ST required to excite one new lattice site. In other words, ST(L) is the time required 
to pass from L to L + 1 lattice sites, or the first passage time of the spreading. We define 
a lattice site as excited when its local energy exceeds some border, in numerics below we 
take Eb = 10~ 50 . The actual value of Eg was chosen arbitrarily, but any other value, 
e.g. Eb = 1CT 100 would produce similar results. The quantity ST can be interpreted as 
a propagation time for L — > L + 1, it can be determined for each particular realization 
of disorder and initial conditions. After having the ensemble of these propagation times 
at given L, we calculate the average propagation time AT as the geometric average of 
ST (equivalently, we average the logarithms: AT = exp[(log ST)}). 

This second approach is superior to the measurement of the effective spatial extent 
C(t) for two main reasons: 

(i) First, AT being the time interval between two spreading events, has no explicit time 
dependence, hence any prehistory does not appear in later measurements of AT. It is 
therefore easier to compare different realizations and simulations for different parameter 
values and different initial conditions using AT. 

(ii) The second advantage relates to the procedure of averaging over many realizations 
of spreading trajectories. By averaging AT(L) for fixed L, we average over situations 
with the same energy density w = S/L. If, in contrast, different realizations of C(t) 
are averaged for a fixed time t, situations with different energy densities w are averaged 
together, which is not reasonable if the density w is the crucial parameter on which 
the properties of the propagation should depend. The difference between averaging 
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Figure 4. Left: Self similar solution p(x) as given by eq. ((9|) of the NDE for a = 2 at 
times t = 10 4 , 10 5 , 10 6 , 10 7 , 10 s (inner to outer lines). Note the logarithmic scaling of 
p. Right: Averaged local energy excitation (E n ) of the spreading state in a nonrandom 
lattice at times T = 10 6 and T = 10 7 for n — A = 6 and f3 = 4. The average is taken 
over M = 48 random initial conditions. The dashed black line shows the corresponding 
analytic self-similar solution of the NDE. 



procedures for C{t) and AT(L) is illustrated in Fig. |3j 
3. Nonlinear Diffusion Equation 

We study spreading which is induced by nonlinear chaotic interactions between 
oscillators, it disappears in the linear limit. Phenomenologically, we have seen in Fig. [2] 
that the spreading looks like a slow diffusion. While we are not able to construct 
a consistent theory for this process, we will stick to a heuristic approach and use 
a phenomenologic model to describe the spreading process: the Nonlinear Diffusion 
Equation (NDE): 

dt dx if dx ) a + 1 dx 2 ^ ' I ^ ^ 

Here p(x, t) is the energy density as a function of position x and time t, and a is a 
parameter of the nonlinear exponent, which later will be related to the nonlinear ities 
in the lattice models. The second equation represents the energy conservation in this 
process. 

Our main method is to compare the (average) spreading of energy in our lattice 
model with predictions by the NDE. The NDE has already been successfully applied 
to nonlinear spreading in the DANSE model [14J and Hamiltonian chains like the ones 
discussed here [T6] . 

3.1. Self- similar solution 
The self-similar solutions of NDE reads [22] 





.2 \ l/ a 

~ 2(a+2)(t-t ) 2 A 2 +«) J for \A< X 

for \x\ > X 
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with c being a constant of integration which follows from the energy conservation: 

- ( I a r(3/2 + l/q) \ 2a/(a+2) 
°~ \y2ir(a + 2) ' T(l + 1/a) J 

X is the sharp front of the field and has the following time dependence: 



X = X 2c^^£ a ' {2 + a \t - t ) 1/{2+a) . (10) 
V a 

The solution has sharp edges (see Fig. H|) and its spatial extension is given by X. 
The size of the wave packet grows in time as a power law, which can be represented in 
a scaling form as 

In order to get rid of the undetermined time offset to, we calculate, following [16] . the 
local inverse velocity of the spreading as dt/dX and obtain for it the following scaling 
law: 

ldt (12) 



'11} 



XdX \X 
3.2. Implications for spreading in lattices 

In our numerical simulations of strongly nonlinear lattices, one approach is to calculate 
the characteristic size of the wave packet by appropriate averaging of the entropies 
(J7J). We can directly associate the size of the field £, with the width of the self-similar 
solution X, so that in (JTTT) X ~ C. Thus, if we assume that the NDE provides an 
adequate description of the spreading in strongly nonlinear lattices, the spreading data 
for different energies should fulfill scaling (II ip. where the constant a depends generally 
on the powers k, A: 

Similarly, in the second method we calculate the average time AT needed for 
spreading, in dependence of the field spatial extend L (which is equivalent to the size 
X of the NDE self-similar solution). This quantity directly corresponds to the inverse 
velocity in (IT21 : dt/dX ~ AT. Thus, we expect that the times AT scale as 

AT /£y a{w) 

Here we introduced the density w as w = S/L, and have generalized the scaling, allowing 
for a possible (slow) dependence of the constant a on this density (what will be confirmed 
by numerics below). 



(14) 
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3.3. Scaling for homogeneously nonlinear lattices 

In the case of lattices (CQ) with k = A, i.e. when both on-site and coupling nonlinearities 
have the same power, we can find the relation between the index a and parameter k 
analytically. Indeed, from the scaling (j3J) we obtain that the time scales with energy as 

t~£^. (15) 

Comparing this with the scaling t — to ~ X 2+a £~ a we obtain 

a = —— , X~(t-t )**, — ~X 1+a . (16) 

IK CJA 

These exact relations will serve as a test for our assumption that the NDE (JH} adequately 
describes the spreading in lattices (TTT) . 



4. Numerical Results 



In the following sections we report on extensive numerical simulations of strongly 
nonlinear lattices, trying to check the predictions of the NDE framework. For the 
numerical time evolution we used a 6th-order symplectic Runge-Kutta scheme [23J, 
mostly with step-size At = 0.1. 



4-1. Homogeneous nonlinearity 

We start with the case of homogeneous nonlinearity k — A, where the scaling index 
follows from relations ({TBI . The Hamiltonian flH]) contains one dimensionless parameter 
(3, additionally the distribution of the local frequencies Uk should be defined. To test 
the theoretical predictions ( fl6l) we follow the evolution in a one-dimensional lattice with 
u>k G [0, 1], started from a single site (or several sites for k = 6) excitation in the middle. 
For several values for the nonlinear strength /3 = 1/4,1/2,1,2 we we integrated the 
system up to T = 10 6 and analyzed the spreading in terms of P(t) and AT(L). This 
was repeated for M = 100 realizations of disorder. Fig. [5] shows the averaged results of 
these runs for the excitation time AT(L) for k = 4 (left panel) and k = 6 (right panel). 
In both cases we find a quite nice agreement of the numerical results with the analytic 
predictions of the NDE (fTBT) . We also performed simulations choosing the disorder to 
be Uk G [0.5, 1.5] and got similar results (not presented here). Additionally, results for 
the direct spreading measure P(t) show the same agreement with prediction ffl6|) (see 
ref. [16J). To our opinion, the agreement between numerics and prediction is a rather 
convincing evidence that the NDE provides the proper framework to model the average 
energy spreading in nonlinear Hamiltonian chains. 

Remarkably, the NDE scaling ( fl6l) holds also for one-dimensional lattices without 
disorder. Fig. [6] shows the participation number evolution for lattice (jSJ) with Uk = 1 
and k = 4,6. In this case the excitation times AT are not the proper measures as 
they are dominated by quasi-compactons (cf. Fig. [D(a)), but the participation number 
calculations are insensitive to such modes. For them we expect from (Tl6l) the scaling 
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Figure 5. Excitation times AT for the lattices with homogeneous nonlinearities, 
k = A = 4 (left panels) and k = A = 6 (right panels). The black lines show the 
predicted behaviors [see (TIB)) ] AT ~ L 5 / 4 (left) and AT ~ L 4 / 3 (right) respectively. 



P ~ [t — t ) 5 «- 2 which is confirmed in Fig. [6j This is also supported by the direct 
comparison of spreading states from the numerical simulation to the self similar solution, 
as can be seen in right panel of Fig. HI 




Figure 6. Participation number for a nonrandom lattice with k = A = 4 (left panel) 
and k = A = 6 (right panel) averaged over M = 48 random initial conditions (solid 
lines). The dashed lines shows the NDE prediction P ~ i 4 / 9 and P ~ t 3 / 7 respectively. 
In the inset we plot the numerical spreading exponent v obtained from finite differences 
of the method above, the dashed line there also corresponds to the expectation from 
the NDE. Notice that a slight final deviation of the spreading exponent for j3 = 1, 2 in 
the inset of Fig. [6] is due to boundary effects as we hit the maximum lattice size at the 
end of integration. 



4-2. Linear on-site Potential, Nonlinear Coupling 

After having found that the NDE provides a good description of the spreading for the 
case of homogeneous nonlinearity, we turn now to a general situation with k ^ A. In 
this section we assume the on-site potential to be quadratic k = 2, for the coupling we 
chose A = 4 and A = 6. This is the case of nonlinear ly coupled linear oscillators. It 
corresponds to a rather general situation of nonlinear disordered lattices, where in the 
representation of linear eigenmodes one can also interprete the system as an ensemble of 
nonlinearly coupled linear modes. The main difference of our setup of strongly nonlinear 
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Figure 7. Excitation times AT(L) for the case with linear on-site term and nonlinear 
coupling k — 2, A = 4. The upper panels show results for disorder chosen as LOk = [0, 1] 
while in the lower two graphs we had ujk = [0.5, 1.5]. In the left panels you see direct 
results AT(L) while on the right the scaling from the NDE has been applied and we 
plot AT/L vs. L/£. Each color/symbol belongs to an averaged value over disorder 
realization for a fixed energy £ . The inlets on right panels show dependence of the 
nonlinearity index a(w) on the density w, obtained via polynomial fitting of the data 
(dashed red lines). 



lattices from that of coupled linear modes, is that here the linear oscillators interact 
only with the nearest neighbors, and the interaction terms contain no disorder, while in 
general lattices, the modes interact with several other modes and the interaction terms 
are disordered, because the coupling involves the overlap of the exponentially localized, 
disordered eigenmodes. 

As outlined in section 12.21 the only free parameter in the Hamiltonian ((3]) is the 
conserved energy S, and we analyze the excitation times AT(L) as function of the wave 
packet length L, for different energies S. Our goal is to check the predictions of the 
NDE scaling (jlijl . 

In Fig. [7] we report the results for k = 2 and A = 4. We have performed many 
runs for different realizations of disorder choosing ujk G [0, 1] and Uk G [0.5, 1.5], and for 
different total energies (some of these results have already been presented in [16]). For 
each disorder realization, several sites are initially excited, and from the evolution of 



Scaling of Energy Spreading in Strongly Nonlinear Lattices 



12 



the field the propagation times AT(L) were determined as described in section [231 Left 
panels of Fig. [7J show these times in dependence on L for different energies. Applying 
the energy scaling according to ( fill) , however, we find a perfect collapse to a single curve 
for all energy values, as demonstrated in right panels of Fig. [7J The same collapse of 
data can be observed in Fig. [HI for k = 2, A = 6, where we only show the rescaled plots. 

Hence, we found that for both cases A = 4, 6 and the two disorder ranges [0,1], 
[0.5, 1.5], the energy scaling resulting from the NDE is valid over more than two orders 
of magnitude of the total energy £. However, from the scaled plots one can see that 
the resulting graph is not a straight line with slope a as predicted by the NDE, but a 
curved one. This means that the nonlinearity index a is effectively density dependent, 
so we define 

d (AT\ 
a{w) = ~fa {—) 

and show this values in the inlets. The curves bend upwards and hence a(w) increases 
with decreasing energy density (note that the x-axis in Figs. [7J and M represent L/S = 
w^ 1 ) resulting in longer excitation times and slower spreading for smaller densities. 

We conclude that in the 2, A = 4, 6 we find the NDE to give an 

adequate description in the sense of correct energy scaling. However, numerical results 
here indicate that the nonlinearity parameter in the NDE is (logarithmically) density 
dependent a(w), and increases with decreasing density. As the density decreases in 
course of the excitation spreading, this also results in a decrease of the effective spreading 
exponent (2+a) -1 with time. The question of the asymptotic behavior remains, however, 
open: from the data presented we cannot judge whether the spreading effectively stops, 
or continues with an increasing index a, or some transition to another law of spreading 
(e.g., a logarithmic one) occurs. 
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Figure 8. Scaled excitation times AT(L) for the lattice with linear on-site term and 
nonlinear coupling k — 2, A = 6. In the left panels you see the scaled results AT/ L 
vs. L/£ for u!k S [0,1] while on the right we set uik & [0.5,1.5]. Each color/symbol 
belongs to an averaged value over disorder realization for a fixed energy. The inlets 
show the density dependencies a(w) obtained through the cubic approximation of the 
data points (dashed red curves). 
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Figure 9. Excitation times AT(L) for k — 4, A = 6 with energies £ = 0.001 . . .0.05. 
The left panel shows the original data while on the right panel scaled variables as 
suggested by the NDE are used: AT/L vs. L/£ . The dashed black line has slope 1.66. 



4-3. Nonlinear on-site Potential, Nonlinear Coupling 

In this section we describe the lattice with A = 4 and k = 6, i.e. a lattice of nonlinearly 
coupled nonlinear on-site oscillators. In this case, the disorder parameter Uk in ([3]) does 
not have the meaning of an oscillator frequency, but is the coefficient determining the 
nonlinearity strength. The real frequency of oscillations depends on the local energy at 
the site. In Fig. [T] we show an exemplary time evolution of an initially localized state in 
such a lattice. 

Again, we investigate the excitation times AT(L) for different total energies in the 
same way as above. Results of our simulations for ujf. G [0.5, 1.5] are shown in Fig. [9j 
In the right panel the scaling as suggested by the NDE ([I4~]) is applied. We again 
find a perfect overlapping of the numerical curves supporting the validity of the NDE. 
Contrary to the previous result for nonlinear coupled linear oscillators (« = 2), we find 
that data of the right panel of Fig. [9] can be nicely fitted with a straight line. This means 
that the parameter a does not dependent on the density. Hence, in this case the NDE 
gives not only the correct scaling, but can be also applied in the standard form with a 
constant effective parameter a. The results in Fig. [9] show roughly AT/L ~ (L/S) im 
which means, according to ffl~2l . a ~ |. This gives for the participation number the 
subdiffusive spreading P ~ t 3 ^ 11 . In Fig. [10] we check this by directly analyzing P{t) 
for different energies. As one can see in the right panel, the scaling predicted from the 
NDE (TiTf works rather well for an enormous range of energy values £ = 0.001 ... 1. The 
scaling is not as perfect as for the excitation times above, which we believe comes from 
the fact that the excitation times are better representatives of the spreading behavior in 
this system as discussed in section 12.41 (and, additionally, no fitting of the time offset to 
has been performed). The asymptotic slope of the scaled curve is in a good accordance 
with the results for the excitation times in Fig. [91 
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Figure 10. Participation Number P(t) for k — 4, A = 6 and energies £ = 0.001 . . . 1. 
The left panel shows original data while on the right the scaling variables P/£ vs. t/£ 2 
are shown. The slope of the red line is 3/11 which corresponds to a = 5/3 as found in 
Fig. [9] above. 



5. Conclusions 

Motivated by previous observations of subdiffusive behavior in nonlinear disordered 
systems, we used the Nonlinear Diffusion Equation as a phenomenologic model to 
describe the spreading process in disordered one-dimensional Hamiltonian lattices of 
nonlinearly coupled oscillators. We have found that the NDE not only predicts a 
subdiffusive spreading, but also induces a one-parameter scaling of time and energy 
of the spreading process according to relations ( !T3lfT4l) . We tested these scaling laws on 
a class of nonlinearly coupled oscillators with different values of the nonlinear exponents 
k and A. In all considered cases the one-parameter scaling laws were nicely reproduced in 
our numerical study. In this sense one can say that the spreading in nonlinear disordered 
one- dimensional lattices belongs to a "universality class" of NDE. For nonlinearly 
coupled linear disordered oscillators, we have found a remarkable dependence of the 
effective index of nonlinearity of NDE on the energy density. As the effective nonlinearity 
increases in the course of spreading, this means a slowing down of the spreading process 
compared with the perfect power law. Still, it is not possible to judge from the data, 
what is the asymptotic regime of the spreading. 

The essential step of our study is an introduction of a novel quantity, characterizing 
the spreading - the averaged excitation time AT. This quantity, defined as a first 
passage time for a particular size of the wave packet, and thus for a particular value of 
the density, allows us to reveal the density dependence of the spreading characteristics, 
what is hardly possible with old approaches where, e.g., the averaged participation 
numbers have been followed (cf. a clear scaling in Fig. |9] with a weakly converging 
asymptotic in Fig. [TO]) . Unfortunately, the calculation of the averaged excitation time 
is relied on the sharp edges of the field, so its application to linearly coupled lattices 
where eigenmodes are exponentially (but not sharp) localized, remains a challenge for 
future studies. Also, in future studies it would be desirable to extend the analysis to 
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two-dimensional lattices. 
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